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ABSTRACT 



We describe the processing of data from the Low Frequency Instrument (LFI) used in production of the Planck Early Release 
Compact Source Catalogue (ERCSC). In particular, we discuss the steps involved in reducing the data from telemetry packets to 
k^. , cleaned, calibrated, time-ordered data (TOD) and frequency maps. Data are continuously calibrated using the modulation of the 
f— ^ • temperature of the cosmic microwave background radiation induced by the motion of the spacecraft. Noise properties are estimated 
from TOD from which the sky signal has been removed using a generalized least square map-making algorithm. Measured 1// 
noise knee-frequencies range from ~ 100 mHz at 30 GHz to a few tens of mHz at 70 GHz. A destriping code (Madam) is employed to 
combine radiometric data and pointing information into sky maps, minimizing the variance of correlated noise. Noise covariance 
matrices required to compute statistical uncertainties on LFI and Planck products are also produced. Main beams are estimated 
down to the « — 10 dB level using Jupiter transits, which are also used for geometrical calibration of the focal plane. 

Key words, methods: data analysis - cosmic background radiation - cosmology: observations - surveys 



.,— I \ 1. Introduction High Frequency Instrument (HFI) (Lamarre et al. 2010; 

■ Planck HFI Core Team 2011a) covers the 100, 143, 217, 

u ■ Planck (Tauber etal. 2010; Planck Collaboration 2011a) is 353j 545j and 857 GHz bands with bolometers cooled to 

_C3 ; a third generation space mission to measure the anisotropy .1 K. Polarization is measured in all but the highest two 

of the Cosmic Microwave Background (CMB). It observes bands (Leahy et al. 2010; Rosset et al. 2010). A combi- 

the sky in nine frequency bands covering 30-857 GHz with nation of ra di a tive cooling and three mechanical coolers 

high sensitivity and angular resolution from 31'to 5'. The pro duces the temperatures needed for the detectors and 

Low Frequency Instrument (LFI) (Mandolesi et al. 2010; optics (pi anc k Collaboration 2011b). Two Data Processing 

Bersanelli et al. 2010; Mennella et al. 2011) covers the 30, Centres (DPCs), conceived as interacting and complemen- 

44, and 70 GHz bands with amplifiers cooled to 20 K. The tary since the earliest design of the Planck scientific ground 

segment (Pasian & Gispcrt 2000); check and calibrate the 



1 Planck (http://www.esa.int/PZancfc) is a project of the data and make maps of the sky, this paper and (Planck 

European Space Agency (ESA) with instruments provided by HFI Core Team 2011b). Planck's sensitivity, angular reso- 

two scientific consortia funded by ESA member states (in par- lution, and frequency coverage make it a powerful instru- 

ticular the lead countries France and Italy), with contributions me nt for galactic and extragalactic astrophysics, as well as 

from NASA (USA) and telescope reflectors provided by a collab- cosmology. Early astrophysics results are given in Planck 

oration between ESA and a scientific consortium led and funded Collaboration 2011c-z 
by Denmark. 
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The Low Frequency Instrument LFI on Planck com- 
prises a set of 11 radiometer chain assemblies (RCAs), each 
composed of two independent, pseudo-correlation radiome- 
ters. There arc two RCAs at 30 GHz, three at 44 GHz, 
and six at 70 GHz. Each radiometer has two independent 
diodes for detection, integration, and conversion from ra- 
dio frequency signal to DC voltage. The LFI is cryogenically 
cooled to 20 K to reduce noise, while the pseudo-correlation 
design with reference loads at w 4K ensures good suppres- 
sion of 1// noise (Mcnnclla ct al. 2011). 

LFI produces full-sky maps centered near 30, 44, and 
70 GHz, with significant improvements with respect to cur- 
rent CMB data in the same frequency range. Careful data 
processing is required in order to realize the full potential 
of LFI and the ambitious science goals of Planck, which 
require that systematic effects be limited to a few /iK per 
resolution element. 

In this paper we describe the processing steps imple- 
mented to create LFI data products, with particular atten- 
tion to the needs of the first set of astrophysics results. 

The structure of the paper follows the flow of the data 
through the analysis pipeline. Section 2 describes the cre- 
ation of time ordered information (TOI) from telemetry 
packets, time stamping, pointing reconstruction, and data 
flagging. Section 3 describes the main operations performed 
on the TOI, including removal of frequency spikes, creation 
of differenced data, determination of the gain modulation 
factor, and diode combination. Beam reconstruction is dis- 
cussed in Sect. 4, calibration in Sect. 5, and noise in Sect. 6. 
Map-making, covariancc matrices, and tests based on jack- 
knife analysis and Monte Carlo simulations arc described in 
Sect. 7. Section 8 reports on colour corrections. Section 9 
describes how the CMB was removed from LFI and HFI 
maps. Finally, Sect. 10 gives an overview of the software 
infrastructure at the LFI DPC. 



2. Creation of time ordered information 

The task of the Level 1 DPC pipeline is to retrieve all neces- 
sary information from packets received each day from the 
Mission Operations Centre (MOC) and to transform the 
scientific TOI and housekeeping (H/K) data into a form 
that is manageable by the scientific pipeline. 

During the ~3h daily telecommunication period 
(DTCP), the MOC receives telemetry from the previous 
day, archived in on-board mass memory, together with real- 
time telemetry. Additional auxiliary files, such as the atti- 
tude history file (AHF) of the satellite, are produced. 

The MOC consolidates the data for each day, checking 
for gaps or corrupted telemetry packets, then provides the 
data, together with additional auxiliary data, to the DPCs 
through a client /server application called the data disposi- 
tion system (DDS). 

The data are received at the DPC as a stream of packets, 
which arc handled automatically by four Level 1 pipelines: 
Data Receipt, Telemetry Handling, Auxiliary Data, and 
Command History. 

The Data Receipt pipeline implements the client side of 
the interface with the DDS. It requests a subset of data pro- 
vided through this interface. A finite-state machine model 
has been used in the design of this pipeline for better for- 
malization of the actions required during interaction with 
the DDS server. 



The Telemetry Handling pipeline is triggered when a 
new segment of telemetry data is received. The first task 
(Telemetry Unscramblcr) discriminates between scientific 
and housekeeping telemetry packets. Scientific packets are 
grouped according to radiometer, detector source, and pro- 
cessing type, then uncompressed and decoded (see next 
paragraph). The on-board time of each sample is computed 
based on the packet on-board time and the detector sam- 
pling frequency. Housekeeping telemetry packets are also 
grouped according to packet type, and each housekeeping 
parameter within the packet is extracted and saved into 
TOI. Subsequent tasks of the pipeline perform calibration 
of housekeeping and scientific TOIs together with addi- 
tional quality checks (e.g., out of limits, time correlation). 
The last task, FITS2DMC, ingests the TOIs into the Data 
Management Component (DMC), making them available 
to the Level 2 and Level 3 pipelines. 

The Auxiliary Data pipeline ingests the AHF provided 
by Flight Dynamics into the DMC. Finally the Command 
History pipeline requests and stores the list of telecom- 
mands sent to the satellite during the DTCP. 

The four pipelines are implemented as perl scripts, 
scheduled every 5 min. Trigger files are created to activate 
the processing in the Auxiliary Data and Command History 
pipelines, and a pipeline monitoring facility displays infor- 
mation about the status of each pipeline. The entire Level 1 
pipeline was heavily tested and validated before the start of 
Planck operations (see Frailis et al. 2009, for more details). 



2.1. Scientific data processing 

When creating TOI, the Level 1 pipeline must recover ac- 
curately the values of the original (averaged) sky and load 
samples acquired on-board. The instrument can acquire sci- 
entific data in several modes or "PTypes" ; we describe here 
only the nominal one (PType 5) (see Zacchei et al. 2009). 
The key feature is that two independent differenced time 
streams are created from the sky and load signals with two 
different gain modulation factors (GMFs). 

Data of PType 5 are first uncompressed. The lossless 
compression applied on-board is inverted, and the number 
of samples obtained is checked against auxiliary packet in- 
formation. Decompressed data Qj-1,2 are then subject to 
a dequantization step to recover the original signals Pi ac- 
cording to 



P 



Qi 



SECOND.QUANT 



OFFSET_ADJUST , (1) 



where SECOND.QUANT and OFFSET_ADJUST are pa- 
rameters of the readout electronics box assembly (REBA), 
calibration of which is described by Maris et al. (2009). 

After dequantization, data are demixed to obtain S^ky 
and S'load using as inputs the gain modulation factors R\ 
and i?2 determined during REBA calibration (Maris et al. 
2009): 

- R 2 ■ P 1 - R 1 ■ P 2 

b sky - 5 5 , {<!) 
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Pi-P* 

i?2 — Rl 



(3) 



2 



Zacchei A., et al.: LFI data processing 



Conversion from ADU (analog-to-digital units) to volts 
is achieved by 



Vi = Sl ~ ZDAE Odae ■ 

LrDAE 



(4) 



where Gdae, Odae, and Zdae are data acquisition elec- 
tronics (DAE) gain, offset, and small tunable offset, re- 
spectively, whose optimal values were determined during 
ground tests (Maris et al. 2009). 

2.2. On-Board Time reconstruction 

A time stamp is assigned to each data sample. If the phase 
switch (Mennella et al. 2011) is off (not switching), the 
packet contains consecutive values of either sky or load sam- 
ples. Then 



t° 



bt 



'() 



bt 



f 



(5) 



where 

xobt 



snip 



> is the sample index within the packet and 
tg bt is the mean time stamp of the first averaged sample. 
•Waver is the number of fast samples averaged together to 
obtain a single detector sample, and / samp ~ 4 kHz is the 
detector sampling frequency. 

If the phase switch is on (nominal case), consecutive 
pairs of cither sky— load or load— sky samples are stored 
in the packet. Then consecutive pairs of samples have the 
same time stamp and 



bt 



bt 



2 trunc 



^smp 



(6) 



On-board time information is stored in the form of TOI 
and directly linked to its scientific sample. 

2.3. Data flagging 

For each sample we define a 32-bit flag mask to identify po- 
tential inconsistencies in the data and to enable the pipeline 
to skip that sample or handle it differently during further 
processing. Currently flags that are checked include: those 
identifying the stable pointing period (determined from the 
AHF); science data that cannot be recovered (e.g., because 
of saturation); samples artificially created to fill data gaps; 
and samples affected by planet transits and moving objects 
within the Solar system. 



3. TOI processing 

3.1. Electronic spikes 

The clock in the housekeeping electronics is inadequately 
shielded from the data lines, resulting in noise "spikes" in 
the frequency domain at multiples of 1 Hz (Meinhold et al. 
2009). The spikes are synchronous with the on-board time, 
with no change in phase over the entire survey, allowing 
construction of a piecewisc-continuous template by sum- 
ming the data for a given detector onto a one second interval 
(Fig. 1) The amplitude and shape differ from detector to de- 
tector; differences between detectors of different frequency 
tend to be larger than between those of the same frequency. 
The amplitude also varies with time. This variation is esti- 
mated by constructing templates like Fig. 1 summed over 
the entire survey to obtain the shape of the signal, and then 



fitting the amplitude of a signal of this shape for each hour 
of data. This amplitude is smoothed with a 20-day boxcar 
window function to reduce the noise. Because of noise, this 
is likely to be an overestimate of the true variations. 
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Fig. 1. A square wave template for both sky (black) and 
load (red/dashed line) for one of the 44 GHz detectors, com- 
puted by adding data between Operational Day (OD) 91 
and 389 in phase over a 1-hr interval. Individual templates 
are directly subtracted from the un-diffcrcnccd data. 



To estimate the effect of spikes on the science data, we 
generate three simulated maps at each frequency. The first 
is a noise map, generated from the instrument white noise 
levels as measured in the data and the scan strategy of 
Planck, but no spikes or correlated noise. This is a best 
case scenario, with the lowest noise level possible in a real 
map. The second is a "spike" map, calculated assuming 
the square wave template for each detector modulated by 
a time- varying amplitude measured from the data, as de- 
scribed above. Because the variation of amplitudes is an 
overestimate, as described above, this is a worst-case sce- 
nario of the effect of spikes. The third map is a "spike- 
subtracted" map, the same as the second, but with a con- 
stant spike template subtracted. This gives an estimate of 
the residual effect of electronic spikes that would be left in 
the maps if the spike template were subtracted. The 30 GHz 
maps are at HEALPix resolution N s ia c = 512; the 44 GHz 
and 70 GHz maps were produced at N s id e = 1024. 

We scale the spike and spike subtracted maps to the 
noise, i.e., Map2/rms(Map3) and Map2/rms(Map3), where 
the rms is calculated from the global rms of the noise map 
scaled as appropriate for the relative number of observa- 
tions ("hits") in that pixel. Figure 2 (left) shows the max- 
imum value of these ratios over the whole sky. At 44 GHz, 
the most affected frequency, the effect in the worst pixel 
is less than 20% of the noise. At 70 GHz the effect in the 
worst pixel is an insignificant 2% of the noise. 

Figure 2 (right) gives angular power spectra of the three 
44 GHz maps. The effect is everywhere well below the noise, 
and subtraction of a constant amplitude square-wave tem- 
plate reduces the effect by almost two orders of magnitude. 

We decided to remove a square-wave template only at 
44 GHz. This reduces the spike residual from 20% of the 
noise to 1% of the noise. At 70 GHz the effect of spikes is 
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Fig. 2. The effect of electronic spikes on the data. Left — Maximum pixel value in the simulated "spike" (black dots) 
and "spike-subtracted" maps (red triangles; see text), scaled to the local pixel noise. In our data processing, the square 
wave signal is subtracted only at 44 GHz detectors. The black circles therefore represent the estimated highest spike 
signal level in the 30 GHz and 70 GHz maps, while the red triangle represents the estimated highest residual spike signal 
level in the 44 GHz map. Right — Angular power spectra of the 44 GHz simulations: the red (middle) line shows the 
power spectrum of the simulated spike map, and the green (bottom) line shows the power spectrum of the simulated 
spike-subtracted map. Subtraction reduces the power by a factor of about 100, from a small to an insignificant fraction 
of the white noise power, shown by the black (top) line. 



extremely small without correction, and at 30 GHz uncer- 
tainty in the template combined with the small size of the 
effect argued against removal. 

3.2. Gain modulation factor and differenced data 

The output of each detector (diode) switches at 4096 Hz 
(Mennella et al. 2010) between the sky, V^ky, and the 4K 
reference load, Vioad- Kky and Vi oa d are dominated by 1// 
noise, with knee frequencies of tens of hertz. This noise 
is highly correlated between the two streams, a result of 
the pseudo-correlation design (Bersanelli et al. 2010), and 
differencing the streams results in a dramatic reduction of 
the 1// noise. The two arms of the radiometer are slightly 
unbalanced, as one looks at the 2.7 K sky and the other 
looks at the ~ 4.5 K reference load. To force the mean of 
the difference to zero, the load signal is multiplied by the 
GMF, R, which can be computed in several ways (Mennella 
et al. 2003). The simplest method, and the one implemented 
in the processing pipeline, is to take the ratio of DC levels 
from sky and load outputs obtained by averaging the two 
time streams, i.e., R = (14ky)/(Moad)- Then 

AV(t)=V sky (t)-^-V load (t). (7) 

\ Vioad/ 

We compute R from unflagged data for each pointing period 
identified from the AHF information. 

To verify the accuracy of this approach, we started with 
a time stream of real differenced data, then generated two 
time streams of undifferenced data using a constant (typi- 
cal) value of R. We then ran these two time streams through 
the pipeline, and compared the results with the original 
time stream. Deviations between the pipeline values of R 
and the constant input value used to generate the undiffer- 
enced data were at the 0.002% level. 



The R factor has been stable over the mission so far, 
with overall variations of 0.03-0.04%. To keep the pipeline 
simple, we apply a single value of R to each pointing period. 

Figure 3 shows the effect of applying Eq.7 with the R 
factor to flight data. The correlated 1// noise in sky and 
load streams (evident in the two upper plots of the figure) 
is reduced dramatically. The residual 1// noise has a knee 
frequencies of 25 mHz, and little effect on maps of the sky, 
as described in Sect. 7. 

3.3. The diode combination 

Having two diodes for each radiometer enables observation 
of both sky and load with a combined duty cycle of almost 
100%. In combining the outputs, however, we must take 
into account the effects of imperfect isolation and differ- 
ences in noise between the two diodes. 

Isolation between diodes was measured for each ra- 
diometer in ground tests and verified in flight using the 
CMB dipolc, planets, and Galactic plane crossings. Typical 
values range from —13 to — 20 dB. This is within specifi- 
cations, and does not compromise LFI sensitivity. It does, 
however, produce a small anti-correlation of the white noise 
of the two diodes of a given radiometer. When data from 
the two diodes are averaged, the white noise of the resulting 
TOI is lower than would be the case if they were statistically 
independent. A complete mathematical description of this 
behaviour is given in Mennella et al. (2011). This causes no 
difficulty in subsequent calibration and further processing; 
however, the effect must be taken into account in inferring 
the noise properties of individual detector chains from the 
combined outputs. 

To take account of differences in noise in combining 
the diode outputs, we assign relative weights to the uncali- 
brated diode time-streams based on their calibrated noise. 
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Fig. 3. Effect of the gain modulation factor (GMF) on sky and load signals for flight data. The upper and middle 
panels show lOmin of sky and load signals of the LFI27S-11 detector: they are highly correlated with clear signatures 
of low- frequency noise. After application of the GMF in taking the difference (Eq. 7), such fluctuations are dramatically 
reduced, revealing the presence of a sky signal dominated by the CMB dipole (lower panel). Note the change in the y-axis 
scale. 



Specifically, we make a first order calibration of the time- 
lines, Go and G\, subtract a signal estimate, and calculate 
the calibrated white noise levels, Co and <ti, for the two 
diodes. The weights for the two diodes (i = or 1) are 



TF, 



1 



G i CTq + a \ 



(8) 



where the weighted calibration constant Gqi is given by 



1 



[G aj + G^l] 



(9) 



and is the same for each diode. 

The weights are fixed to a single value per radiometer for 
the entire dataset. Since all calibrations, noise estimation, 
and other tests are done on these combined data streams, 
small errors in the weights cause inconsequential losses in 
sensitivity, and no systematic errors. 



times for a single pointing period. It specifies with appro- 
priate flags the periods of spin-axis maneuvers during which 
star tracker positions are unreliable. Horn locations within 
the focal plane are determined from both ground measure- 
ments and planet crossings. 

The orientation of the spacecraft spin axis at the time 
of each data sample is determined by linear spherical inter- 
polation of the 8 Hz quaternions. Individual detector point- 
ings are determined by simple rotations from the spin-axis 
reference frame to the telescope optical axis, then to the rel- 
evant horn position, with an additional rotation to account 
for the orientation of the horn in the focal plane. 

In some cases small extrapolations of the quaternions 
are necessary at the end of a pointing period. Simulations 
verify that these introduce no significant degradation of the 
pointing accuracy. 



3.4. Detector pointing 

Detector pointing is a fundamental ingredient in data pro- 
cessing that requires knowledge of the spacecraft attitude 
and the location of the horns in the focal plane. The AHF 
gives the orientation of the spacecraft spin axis in quater- 
nions sampled at 8 Hz, as well as beginning and ending 



4. Main beams and the geometrical calibration of 
the focal plane 

Knowledge of the beams is of paramount importance in 
CMB experiments. Errors and uncertainties, and the de- 
tails of complex non-Gaussian shapes, directly affect cos- 
mological parameters. 
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We determine the main beam parameters and the po- 
sition of each horn in the focal plane from planet observa- 
tions. Jupiter gives the best results, but other planets and 
bright celestial sources have been used as well. Inputs to the 
calculations include TOI from each radiometer throughout 
the planet crossing, the AHF for the same period, and the 
time-dependent position of the planet as seen by Planck, 
provided by the JPL Horizons system, which accounts for 
both spacecraft and planet motion. 

4.1. Algorithm and testing 

We create a 2D map of the footprint of the focal plane on 
the sky by selecting data within 10° of the telescope line 
of sight. This comprises the whole extension of the LFI 
focal plane. To minimize the effects of 1/f noise on weak 
sources, we use TOI from which offsets per ring derived by 
the Madam destriper (Sect. 7) have already been removed. 

We fit a bivariate Gaussian beam model to these data 
(Burigana et al. 2001): 



B{xi,yi) 



exp 



(Axi cosa + Ayi sina) 2 



(—Axi sina + Ay; cosa) 2 



(10) 



Here A is an overall amplitude. Xi and y^ are Cartesian 
coordinates, with x c , y c the position of the centre of the 
beam, and Axi = Xi — x c and Ay; = y^ — y c . o x and <j v are 
the beamwidth parameters of the elliptical approximation 
of the beam shape, and the angle a is the reconstructed 
orientation of the beam in the focal plane and d is the 
actual distance (in astronomical units) of the planet. 

We tested our technique with simulations using the mea- 
sured beam patterns together with a detailed model of the 
Planck telescope. The simulations included the nominal 
main and far beam patterns, the effects of smearing caused 
by the motion of the satellite, and pointing uncertainties. 
Using these simulations of Jupiter crossings (including in- 
strumental noise and complete sky signal), we are able to 
reconstruct the main beam shape down to — 20 dB and to 
recover the main beam properties at the 1% level or better 
for all LFI beams. Table 1 reports results for the main beam 
properties for a sample of the LFI beams. These figures are 
representative of our expected accuracy for in-flight beam 
and focal plane reconstruction. 

Figure 4 shows the footprint of the LFI focal plane ob- 
tained during the first season of Jupiter observations, from 
24 October to 1 November 2009. Figure 5 shows beam im- 
ages for LFI28M, LFI25M and LFI21M from those observa- 
tions. As expected, all beams are asymmetric but with no 
significant departures from an elliptical shape visible down 
to the ~ — 10 dB level. For lower levels, aberration starts to 
distort the beam response, creating non-elliptical shapes. 

We also constructed a planet mask, including Jupiter, 
Mars, and Saturn, the most luminous planets at LFI fre- 
quencies. The planet mask is radiometer-dependent, since 
each horn observes a planet at different times. The planet 
masking algorithm assigns an appropriate flag to data that 
lie within an ellipse, centred at the position of the planet 
and with an orientation that matches the beam orientation, 
with axes ~ 3 times larger than the beam widths derived 
from beam fitting. These flags are used in the map-making 



Table 1. Simulation of the reconstruction of the beams and 
focal plane geometry from observations of Jupiter, includ- 
ing realistic models of the beams, instrument noise, beam 
smearing, and star tracker uncertainties. 
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Fig. 4. LFI focal plane as determined from the first season 
of Jupiter observations, 24 October to 1 November 2009. 
Contour levels are in dB from the peak. All beams are well 
approximated by an elliptical Gaussian down to the — lOdB 
level. 



and ensuing data analysis to discard samples affected by 
planet transits 



5. Photometric calibration 

5.1. First steps 

The ideal source for photometric calibration, i.e., conver- 
sion of the data from volts to kclvin, should be constant, 
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Fig. 5. Beam map of LFI data around the Jupiter observations (24 October-1 November 2009) for LFI horns LFI28M 
(left), LFI25M (middle) and LFI21M (right). 



perfectly known, present during all observations, and have 
the same frequency spectrum as the CMB. In the frequency 
range of the LFI, the CMB dipole, caused by the motion of 
the Solar system with respect to the CMB reference frame, 
satisfies nearly all of these requirements, lacking only in 
that it is well, but not perfectly, known. The modulation 
induced on the CMB dipole by the orbital motion of Planck 
around the Sun satisfies even this last requirement, and will 
be the ultimate calibration source for the LFI; however, it 
cannot be used effectively until data for a full orbit of the 
Sun are available. For this paper, therefore, we must use 
the CMB dipole. We follow essentially the calibration pro- 
cedure used for the WMAP first year data (Hinshaw et al. 
2003). For the fc th pointing period, the signal from each 
detector can be written as 



AV k = g k {AT sky 



(11) 



where AT s k y is the sky signal, n is the noise, and g k and 
b k are the gain and baseline solution. The dominant sky 
signal on short time scales is the CMB dipole (Galactic 
plane crossings produce a localized spike that is easy to 
exclude). This is modeled as 



AV m (g k ,b k ) = g k (AT d + AT V 



(12) 



where we have considered both the cosmological dipole ATd 
and the modulation from the spacecraft motion AT V . We 
fit for g k and b k for each pointing period k by minimising 



E 



[AV(U) - AV m (U\g k ,b k )] 



nils- 



(13) 



The sum includes unflagged samples within a given pointing 
period k that lie outside a Galactic mask. 

The mask is created from simulations of microwave 
emission provided by the Planck Sky Model (PSM) 2 . Of 
the LFI frequencies, 30 GHz has the strongest diffuse fore- 
ground emission. The mask excludes all pixels that in the 
30 GHz PSM are more than 5 x 10~ 4 times the expected 
mis of the CMB. It also excludes point sources brighter 
than 1 Jy found in a compilation of all radio catalogues 
available at high frequencies (the Planck Input Catalogue, 



The Planck Sky Model is available at: 
http : //www. ape .univ-paris7 . f r/APC_CS/Recherche/Adamis/ 
PSM/psky-en . html 



see Massardi 2006). The Galactic and point source masks 
preserve ~ 82% of the sky. 

The Planck scan strategy is such that the instrument 
field of view describes nearly great circles on the sky. The 
signal mean is therefore almost zero and nearly constant 
from one circle to the next. This reduces the correlation 
between the gain and baseline solutions, a feature also taken 
advantage of by WMAP (Hinshaw et al. 2003). 

As pointed out by Hinshaw et al. (2003) and Cappcllini 
et al. (2003), the largest source of error in Eq. 13 arises 
from unmodelled sky signal AT„ from CMB anisotropy 
and emission from the Galaxy. To correct this, we solve 
iteratively for both g' k and AT^. If g' k is the solution at a 
certain iteration, the next solution is derived using Eq. 13 
with 



AV' = AV - g' k AT' a 



(14) 



where AT' & is the sky signal (minus dipole components) 
estimated from a sky map built from the previous iter- 
ation step. This is repeated to convergence, typically af- 
ter a few tens of iterations. Figure 6 shows the gain er- 
ror induced by unmodeled sky signal in a one-year simula- 
tion of one 30 GHz detector. The simulation includes CMB 
anisotropies, the CMB dipole, and Galactic emission. Gain 
errors in this example are ~ 5% after one iteration. After a 
few tens of iterations, the residual errors are < 0.01% over 
the entire year. 

The algorithm alternates between dipole fitting and 
map-making. Maps are made with (Madam Sect. 7) ignor- 
ing polarisation, with no noise prior and baseline length 
equal to the pointing period length. To improve calibration 
and reduce noise, calibration is performed simultaneously 
for both radiometers of each single horn. In the presence 
of real noise, the actual scatter from one gain solution to 
the other is quite large. Figure 7 shows an example of the 
hourly gain solution (grey line) derived from the iterative 
scheme described above for LFI18M, one of the 70 GHz ra- 
diometers. Apart from the scatter induced by instrument 
noise, the gain solution is quite stable throughout the ob- 
servation period. Around the dipole maxima, typical noise- 
induced variations arc ~ 0.8% (rms). Nonetheless, the sta- 
bility of the gain solution is poor compared to the stability 
of the instrument itself, as indicated by the stability of the 
uncalibrated white noise level of both differenced and un- 
diffcrenced data, or the stability of the total power from 
both sky and load signals. This is particularly evident dur- 
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Fig. 6. Simulation showing convergence of the gain solution for one year of observations of one 30 GHz detector. The 
simulations include CMB dipole(s), CMB anisotropics, and Galactic emission. The input gain was 12.86K/V and the 
sky included all diffuse components as well as nominal instrument noise. The first iteration shows large errors caused by 
Galactic and CMB anisotropy emissions; however, after one hundred iterations, convergence is achieved with an overall 
deviation from the input value of less than 0.01%. The various curves show the solution after 1, 30, 60 and 100 iterations. 



ing the minima of the dipole signal (see Mennella et al. 2011 
for further details). 

There are also specific things that affect the gain solu- 
tion. To the extent that they can be measured and under- 
stood, their effects on the gain can be corrected directly. 
For example, a non-linearity in the analogue-to-digital con- 
verters (ADCs), discovered during data analysis, produces 
a multiplicative effect on the data that is recovered (erro- 
neously) by the calibration pipeline as a gain variation. We 
have developed two independent, complementary methods 
to correct for this. In the first, we calibrate the data using 
the gain solution that follows the induced ADC gain varia- 
tion. In the second, we model the nonlinearity and remove 
the effect at the raw TOI level. 

Alternatively, temperature variations of the amplifiers 
can induce real gain variations on short time scales. For ex- 
ample, during the first 259 days after launch the downlink 
transponder was powered up only for downlinks. This in- 
duced rather sharp daily variations in the temperature and 
gain of the amplifiers in the back-end unit (BEM) . Starting 
on day 259, the transponder has been powered up continu- 
ously, eliminating this source of gain variations. 

In the next section we discuss additional steps taken in 
the calibration procedure to deal with the effects of noise 
and gain changes induced by events such as the transponder 
cycle change. 



5.2. Improving calibration accuracy 

As shown in Fig. 7, the hourly gain solutions are strongly 
affected by noise. To reduce the effects of noise and recover 
more accurately the true and quite stable gains of the in- 
strument, we process the hourly gain solution as follows: 

— calculate running averages of length 5 and 30 days. 
The 5-day averages are still noisy during dipole min- 
ima, while the 30-day averages do not follow real but 
rapid gain changes accurately. 

— further smooth the 5- and 30-day curves with wavelets; 

— use the 30-day wavelet-smoothed curve during dipole 
minima; 

— use the 5-day unsmoothed curve around day 259 (the 
downlink transponder change) to trace real gain varia- 
tions; 

— use the 5-day wavelet smoothed curve elsewhere. 

A typical gain solution is plotted in Fig. 7 as the solid 
black line. From the 5- and 30-day gain curves we infer 
information on the actual gain stability of the instrument as 
the mission progresses, and also on the overall uncertainty 
in the gain reconstruction. Specifically, the rms of the gk 
over a period of N pointings is 



Sg 



Efc=i(gfc - (g))' 

N — 1 



(15) 



where (g) is the average of the N gains. The effect of the 
wavelet smoothing filter is to average over a number of con- 
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Tine [days since 5/14/2009] 

Fig. 7. Hourly gain solution (gray line) from flight data for LFI18M, as derived from our iterative calibration algorithm. 
The gain is quite stable over the observing period, although there is a lot of scatter due to noise, especially during dipole 
minima. The thick black line is the refined gain solution (see text) applied to create calibrated TOI and sky maps. 



secutive pointings. Ignoring the different weights in the av- 
erage, the overall uncertainty can be approximated as 



Table 2. Summary of dipole-based gain statistics 



M 



Efc=i(gfc - (g)Y 

N-l 



(16) 



Table 2 lists the largest statistical uncertainties and 
their associated mean gains out of four time windows 
(days 100-140, 280-320, 205-245, 349-389, the first two cor- 
responding to minimum and the second two to maximum 
dipole response), for the main and side arms of the LFI 
radiometers. In order to provide conservative estimates, we 
have always chosen a value for M corresponding to the 
number of pointings in 5 days, even in cases where a 30- 
day smoothing window was used. Equations 15 and 16 and 
Table 2 are the same as equations 12 and 13 and Table 9 
of Mennella et al. (2011). Peak-to-peak variations in the 
daily gains reach 10% (with mean 7%); however, the rms of 
the smoothed gain solution is generally in the ~ 0.3 — 0.4% 
range. This can be taken as the current level of LFI cali- 
bration accuracy. 

Although the current pipeline provides results ap- 
proaching those expected from the stability of the instru- 
ment, we are working to improve it as much as possible. 
In particular, we would like to trace gain variations on 
time scales shorter than the pointing period. To achieve 
this, we are developing a detailed gain model (currently un- 
der test) based on calibration constants estimated from the 
pipeline and instrument parameters (temperature sensors, 
total power data), see Mennella et al. (2011) for further 
information. 
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6. Noise estimation 

Once data are calibrated, we evaluate the noise proper- 
ties of each radiometer. We select data in chunks of 5 days 
each and then compute noise properties. This is done using 
the roma Iterative Generalized Least Square (IGLS) map- 
making algorithm (Natoli et al. 2001; de Gasperis et al. 
2005) which includes a noise estimation tool based on the 
iterative approach described in Prunet et al. (2001). IGLS 
map-making is time and resource intensive and cannot be 
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run over the whole data set within the current DPC sys- 
tem. However since the TOD length considered here is only 
5 ODs, it is possible to use the roma implementation of 
this algorithm which has a noise estimator built-in. The 
method implemented here is summarized as follows. Model 
the calibrated TOD as 



AT = Pm + n, 



(17) 



where n is the noise vector, and P is a projection matrix 
that relates a map pixel m to a TOD measurement AT. 
We obtain a zeroth order estimate of the signal through a 
rebinned map and then iterate noise and signal estimation: 

fii = AT - Prhi, (18) 

ihi+i = (P T Nr 1 P)- 1 P T Nr 1 AT, (19) 

where Nj is the noise covariance matrix in time domain 
estimated at iteration i. We have verified that convergence 
is reached in a few, usually three, iterations. 

We calculate the Fourier transform of the noise time 
stream (with an FFT algorithm) and fit the resulting spec- 
trum for the three parameters, the white noise level, the 
knee-frequency, and the slope of the 1// noise part: 



p(f) = 4 



(20) 



The white noise level is taken as the average of the last few 
percent of frequency bins. A linear fit to the log-log spec- 
trum low frequency tail gives the slope of the 1// noise. 
The knee frequency, /k , is the frequency at which these two 
straight lines intersect. We tested the accuracy with sim- 
ulations that included sky signal and instrumental noise 
with known properties. The noise properties were recov- 
ered with typical deviations from input values of ~ 10% 
for knee-frequency and slope, and less for white noise level. 
Examples of noise spectra and corresponding fits are shown 
in Fig. 8. 

6.1. Noise constrained realizations and gap filling 

The FFT-based noise power spectrum estimation method 
requires continuity of the noise time stream. As discussed 
in Sect. 3, we identify bad data (e.g., unstable spacecraft 
pointing, data saturation effects) and gaps in the data 
with appropriate flags. We fill in the flagged data with a 
Gaussian noise realization constrained by data outside the 
gap (Hoffman & Ribak 1991). Although in principle this 
method requires a pure noise time stream outside the gap, 
we have verified that given the low signal-to-noise ratio in 
the LFI TOD the procedure is not affected by the signal 
present in the time streams. We fill the gap with Gaussian 
noise whose properties match those of the noise power spec- 
trum computed over the day immediately before the one 
with flagged data. An example is shown in Fig. 9. 

7. The Map-Making pipeline 

7.1. Frequency maps 

The map-making pipeline produces sky maps of tempera- 
ture and polarisation for each frequency channel. It takes as 



input the calibrated timelines and pointing information in 
the form of three angles (6, <fi, ip) describing the orientation 
of the feed horns for each data sample. An essential part of 
the map-making process is the reduction of correlated 1// 
noise, a large part which can be removed by exploiting re- 
dundancies in the scanning strategy. While the underlying 
sky signal remains the same, the observed signal varies due 
to noise. Statistical analysis of the signal variations allows 
one to distinguish between true sky signals and noise. 

Among several map-making codes tested with simulated 
Planck data (see Ashdown et al. 2007a,b, 2009) the LFI 
baseline (Mandolcsi et al. 2010) is to use the Madam de- 
striping code (Maino et al. 2002). The algorithm and the 
underlying theory are described in detail in Keihanen et al. 
(2010); Kurki-Suonio et al. (2009); Keihanen et al. (2005). 
The basic idea is to model the correlated noise component 
by a sequence of constant offsets, called baselines. A key 
parameter in the code is the length of the baseline to be 
fitted to the data. Madam allows the use of an optional noise 
prior, if the noise spectrum can be reliably estimated, which 
further improves the accuracy of the output map. Without 
the noise prior, the optimal baseline length is of the order 
of the satellite spin period (« 1 minute). With an accu- 
rate noise prior, a much shorter baseline can be used. The 
shorter the baseline, the closer the Madam solution will be to 
the optimal Generalized Least Square solution (see Fig. 16 
of Ashdown et al. 2009). 

We are continually improving our knowledge of the in- 
strument and its noise characteristics, and this information 
will eventually be used in the Madam algorithm. However, at 
this stage in the processing we decided to make two simpli- 
fications when running our map-making pipeline: no noise 
prior was used, and all radiometers were weighted equally. 
These choices lead to a simpler and faster map-making al- 
gorithm, which is sufficiently accurate for the Planck Early 
Results and avoids using detailed parameters describing the 
instrument which are under continual revision. 

With these simplifications, the map-making equations 
can be written in a concise form. Technically, we are ne- 
glecting the baseline covariance, C a , and setting the white 
noise variance C n to unity. The basic model behind the 
algorithm is 

AT = Pm + n', (21) 

where AT is the calibrated TOD, P is the pointing matrix, 
m is the pixelized sky map, and n' is the instrumental noise. 
This last term can be written as 



n' = Fb + n . 



(22) 



where b is the vector of unknown base function amplitudes 
and the matrix F projects these amplitudes into the TOD. 
Since Madam uses uniform baselines, the matrix F consists 
of ones and zeros, indicating which TOD sample belongs 
to which baseline. Finally n is a pure white noise stream 
assumed to be statistically independent of the baselines. 

The maximum likelihood solution is obtained by mini- 
mizing 



X 2 = (AT Fb Pm) T (AT - Fb - Pm), 



(23) 



with respect to the quantities b and m. The baseline am- 
plitudes b are detcmincd by solving 



(F T ZF)b = F T Zy, 



(24) 
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Fig. 8. Noise spectra of radiometers LFI18M, LFI22S, LFI25S, and LFI28S) estimated by the noise pipeline (black lines). 
All spectra are well- fit by Eq. 8 with a single knee frequency and slope (red lines) . An excess near 1 mHz is visible in 
LFI25S and LFI28S. This is approximately the bed-switching frequency of the sorption cooler, and the different slopes in 
LFI28S and LFI25S on the low-frequency side of the spectrum are possibly indications of thermal effects on the radiometer 
output. 
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Fig. 9. Gap filling procedure applied to LFI28M for day 239. The upper panel shows the original TOI (black) where a 
step is caused by a DAE gain change that produces saturated data. The (red) lines show the constrained noise realization 
used to replace those data. The lower panel shows a zoom around the position of the step to highlight the consistency of 
the gap filling data with the unflagged part of the TOI. 
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where 



Z = I-P(P T P) _1 P 5 



(25) 



Madam uses an iterative conjugate-gradient method to solve 
Eq. (24). An estimate for the map is finally obtained as 



m 



(P T P)" 1 P T (AT-Fb). 



(26) 



The map m has as many elements as pixels in the sky. 
Each element is a Stokes parameter triplet (I, Q, U) for a 
pixel p. The matrix P T P is a 3x3 block diagonal matrix 
that operates on map space. There is a block for each pixel 
p. A block can only be inverted if the pixel p is sampled with 
a sufficient number of different polarisation directions to 
allow determination of the three Stokes parameters for that 
pixel. This is gauged by the condition number of the block. 
For the present analysis, if the inverse condition number 
rcond (ratio of the smallest to largest eigenvalue) was less 
than 0.01, the pixel p was excluded from the (I, Q, U) map. 

The P T P blocks must be inverted when Eqs. (24) 
and (25) are solved for the baselines. These inversions are 
computed by eigenvalue decomposition. Eigenvalues whose 
magnitudes are less than 10~ 6 times the largest eigenvalue 
are discarded; only the remaining part is inverted. 

For the present analysis, we need only the /-component 
maps at the three LFI nominal frequencies, combining ob- 
servations of all radiometers at a given frequency. Figure 10 
shows (left column) the hit count maps by frequency. In 
addition, we produce maps from horn pairs scanning the 
same path in the sky (see Mennclla et al. 2011, for details 
on the LFI focal plane arrangement). We have produced 
30 GHz maps at HEALPix resolution N sidc = 512, and 44 
and 70 GHz maps at A^idc = 1024. All maps are in the 
NESTED scheme, in Galactic coordinates, with units of 
thermodynamic kclvin. The baseline length in Madam was 



one minute 



7.2. White noise covariance matrices 

If we bin the pure white noise stream n to a map using the 
pointing P, we obtain a binned white noise map, 



w = (P T P)- 1 P T n. 



(27) 



This map is a theoretical concept because we do not have 
access to the radiometer white noise streams. Its covari- 
ance matrix, however, is important because it provides an 
estimate of both the white noise power in each pixel and 
white noise correlations between Stokes parameters at a 
given pixel. 

This white noise covariance matrix (WNC) is computed 
as (Eq. 27) 



ww T ) = (P T P)- 1 (P T C„P)(P T P)" 1 



(28) 



Here C n = (nn T ) , C„ is a matrix that operates in the TOD 
domain, and angle brackets denote the ensemble mean. 
Because the radiometers have independent white noise, C„ 
is diagonal. We assume that each radiometer has a uni- 
form white noise variance cr WN (see Sect. 6), but that each 
radiometer has its own variance. The radiometer own val- 
ues that we used in the WNC computation are reported in 
Mennella et al. (2011). 



3 One minute baselines for 30 GHz, 44 GHz, and 70 GHz are 
1950, 2792, and 4726 samples respectively. 



7.3. Half-ring jackknife noise maps 

For noise estimation purposes we divided the time ordered 
data into two halves and produced jackknife maps as fol- 
lows. Each pointing period lasts typically rj 44min (median 
43.5 min, standard deviation lOmin). Typically, during the 
first 4 min the pointing is unstable, so these data are not 
used for science. During the remaining stable 40 min, each 
horn scans a ring on the sky. This ring consists of scan cir- 
cles. One full scan circle takes 1 min. Therefore, each ring 
has about 40 scan circles. We made half-ring jackknife maps 
ji (and j 2 ) with the same pipeline as described in Sect. 7.1, 
but using stable data only from the first or the second half 
of each pointing period. Specificaly, this is implemented by 
marking the other half of each ring as a gap in the data. 
Madam knows that for any given pointing period the first- 
used scan circle/sample of any half-ring is far apart in time 
(typically 25 min) from the last used scan circle/sample of 
the previous half-ring. 

At each pixel p, the jackknife maps ji and j 2 contain the 
same sky signal (as long no time- varying sources or moving 
objects cross p at the time of observation), since they re- 
sult from the same scanning pattern on the sky. However, 
because of instrumental noise, the maps ji and j 2 are not 
identical. 

We can estimate the sky signal+noise as 



mi +2 (p) = [ii(p)+j a (p)]/2, 



and the noise in map rxii + 2 as 



n 1+a (p) = [ii(p)-ja(p)]/2. 



(29) 



(30) 



This noise map includes noise that is not correlated on 
timescales longer than 20 minutes. In particular, ni + 2 gives 
a good estimate of the white noise in mi + 2- 

However, we are interested in the noise level in the full 
map m, (see Eq. 26). To estimate this, we construct another 
noise map 

„ , \ hip) - hip) , Q1 \ 
n m (p) = p-r — , (31) 



w h it(p) 



with weights 



w M t(p) = 4/hitf u ii(p) 



hiti(p) hit 2 (p) 



(32) 



Here hitf u n(p) = hiti(p)+hit2(p) is the hit count at pixel p 
in the full map m, while hiti and hit2 are the hit counts of 
ji and J2, respectively. The weight factor Whit(p) is equal to 
2 only in those pixels where hiti(p) = hita(p). In a typical 
pixel, hiti(p) will differ slightly from hit2(p) and hence the 
weight factor is Whit(p) > 2. 

Noise maps from half-ring jackknifes are shown in the 
right-hand column of Fig. 10. A detailed comparison of the 
jackknife noise estimates and other noise estimates (WNC, 
noise Monte Carlo; see next section) are presented in the 
LFI instrument paper (Mennella et al. 2011). 

7.4. Noise Monte Carlo simulations 

To check the noise analysis, we produced Monte Carlo 
noise realizations on the "Louhi" supercomputer at "CSC- 
IT Center for Science" in Finland. The simulation takes 
as input estimates of the white noise <twn, knee frequency, 
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Fig. 10. Hit count (left) and noise maps (right) at 30 GHz (top), 44 GHz (middle), and 70 GHz (bottom). The complex 
distribution around the ecliptic poles in the 44 GHz hit map is caused by the location of 44 GHz horns on the focal plane. 
The noise maps are derived from half-ring jackknife tests described in the text. 



and slope of the 1/f noise estimated from the TOD for 
each radiometer (Mennella et al. 2011), as well as satel- 
lite pointing information. Flight pointing was reconstructed 
to machine accuracy using Planck Level-S simulation soft- 
ware Reinecke et al. (2006). For each frequency channel, 
we generated 101 Monte Carlo realizations of the noise, 
simulating white noise and correlated noise (1//) streams 
separately. Maps from these noise streams were produced 
with the map-making pipeline described in Sect. 7.1. For 
each simulated noise map, we computed the corresponding 
binned white noise maps defined in Eq. (27). The produc- 
tion of the binned white noise maps allows us to study the 



residual correlated noise, i.e., the difference between the to- 
tal and binned white noise maps Kurki-Suonio et al. (2009). 
These Monte Carlo simulations were used to test and val- 
idate several approaches to noise estimation described in 
detail in the LFI instrument paper (Mennella et al. 2011). 

8. Colour correction 

The power measured by LFI can be expressed as 

P= § I g{v)&T RJ {v)dv , (33) 
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where G is the overall gain, g{v) is the bandpass, and 
ATrj is the Rayleigh- Jeans brightness temperature sig- 
nal, in the case of LFI calibration procedure, due to the 
CMB dipolc. At a given frequency uq, the overall gain G is 
equal to 2P/ ATrj(v )- For small fluctuations around the 
mean CMB temperature To, the relation between intensity, 
/, Rayleigh- Jeans brightness temperature 7rj, and thermo- 
dynamic temperature T is 



2k B^ Arr , \ ( 9B(v ,T) 



AI(u ) = — ^ATrj^q) = 



dT 



AT. (34) 



The differential black-body spectrum is 



dB(v,T) \ 
dT J 



2k B v 2 



o kv/k B T 



T 



hv/k B T G 

z hv/k B T _ I 



2k B v 2 

5—r]AT{v) ■ 



, (35) 



(36) 



The function 77at(^) is the differential black-body spectrum 
in Rayleigh- Jeans units. With our definition of the overall 
gain G. the bandpasses are normalised such that 



gW)T)t\T{v)<lv = ti at (i/ ) . 



(37) 



Calibration data provide a nominal brightness tempera- 
ture ATrj = (2/G)P; however, this is only exact for 
a monochromatic response. For a non-zero bandwidth, a 
colour correction C(a) is required to convert the brightness 
temperature for emission with a particular spectral index 
a to that of the map: 



C{a)AT RJ {u ) = AT RJ = r] AT {is )AT . 



(38) 



By definition, the colour correction is unity when the source 
observed has a CMB spectrum. Within each LFI band, 
g(v) is well- approximated by a power law with spectral in- 
dex a = 2 - {hv /k B T) 2 /6. 

The general expression for the colour correction is 



C(a) 



Vat(vo) 



f g(v)VAT(v)dy 



g^v/v^dv, (39) 



where we assumed a power-law spectrum with temperature 
spectral index (3 = a — 2. The term in square brackets 
is unity with our normalisation for g{v), but has been in- 
cluded to show that C{a) depends only on the shape and 
not the amplitude of the bandpass. Thus C(a) is indepen- 
dent of G. 

Each detector has a different bandpass, hence its own 
colour correction. We derive approximate colour corrections 
for band-averaged sky maps using bandpasses averaged 
over: (i) the two detectors in each radiometer; (ii) the two 
orthogonally-polarised radiometers behind each feed horn; 
and (iii) the several feed horns in each frequency band. In 
addition, although the bandpass is mainly defined by the 
front-end (Bersanelli et al. 2010), differences between back- 
end bandpasses on a single radiometer are measurable, e.g., 
in the form of /3-dependent residuals in difference images. 

Since the current sky maps have been produced, both 
for pairs of horns and for several horns in each band, 
with calibrated data combined with equal weights, we have 



used an unweighted average of all the contributing band- 
passes for our band-averaged corrections. Using the band- 
pass models given in Zonca et al. (2009) derived from the 
pre-launch calibration campaign, we evaluate the integrals 
in Eq.(39) analytically for several spectral indices. The re- 
sults are given in Table 3. 

At the current stage of the mission and data analy- 
sis, uncertainties in the colour corrections are much smaller 
than those of the gains G; however we aim to reduce the 
calibration error (using the orbital dipole) to below 0.2%. 
Two primary sources of error in C(a) will then need to 
be considered. The first is related to uncertainties in the 
bandpass model (Leahy et al. 2010; Zonca et al. 2009). The 
second arises from the uneven sampling of individual sky 
pixels by the full set of detectors, which causes pixel-to- 
pixcl variations in the colour correction. 

9. CMB removal 

This section was developed in common with HFI (Planck 
HFI Core Team 2011b) and is reported identically in both 
papers. 

In order to facilitate foreground studies with the fre- 
quency maps, a set of maps was constructed with an esti- 
mate of the CMB contribution subtracted from them. The 
steps undertaken in determining that estimate of the CMB 
map, subtracting it from the frequency maps, and charac- 
terising the errors in the subtraction are described below. 

9.1. Masks 

Point source masks were constructed from the source cat- 
alogues produced by the LFI pipeline for each of the LFI 
frequency channel maps. The algorithm used in the pipeline 
to detect the sources was a Mexican-hat wavelet filter. All 
sources detected with a signal-to-noise ratio greater than 5 
were masked with a cut of radius 3 er w 1.27 FWHM of the 
effective beam. A similar process was applied to the HFI 
frequency maps (Planck HFI Core Team 2011b). 

Galactic masks were constructed from the 30 GHz and 
353 GHz frequency channel maps. An estimate of the CMB 
was subtracted from the maps in order not to bias the con- 
struction. The maps were smoothed to a common resolution 
of 5°. The pixels within each mask were chosen to be those 
with values above a threshold value. The threshold values 
were chosen to produce masks with the desired fraction of 
the sky remaining. The point source and Galactic masks 
were provided as additional inputs to the component sepa- 
ration algorithms. 

9.2. Selection of the CMB template 

Six component separation or foreground removal algo- 
rithms were applied to the HFI and LFI frequency channel 
maps to produce CMB maps. They are, in alphabetical or- 
der: 

— AltICA: Internal linear combination (ILC) in the map 
domain; 

— CCA: Bayesian component separation in the map do- 
main; 

— FastMEM: Bayesian component separation in the har- 
monic domain; 

— Nccdlct ILC: ILC in the needlet (wavelet) domain; 
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Table 3. Colour corrections for different input power- law spectral indices 



Spectral Index a 



Detector 


-2.0 




-1.0 


0.0 


1.0 


2.0 


3.0 


4.0 


70 GHz 


















LFI18 


1.054 


1 


028 


1.011 


1.003 


1.003 


1.010 


1.026 


LFI19 


1.170 


1 


113 


1.066 


1.026 


0.994 


0.969 


0.949 


LFI20 


1.122 


1 


079 


1.044 


1.017 


0.997 


0.983 


0.975 


LFI21 


1.087 


1 


053 


1.028 


1.010 


1.000 


0.996 


0.998 


T PT99 


u.y / o 


u 


Q71 

y / 1 


u.y / o 


u.yoo 


1 P1H7 
1 .UU ( 


l.Uoo 


1 .uoo 


LFI23 


1.015 


1 


004 


0.999 


0.998 


1.003 


1.012 


1.026 


(C)m 

\W'U 


1 070 


1 


041 


1.021 


1.007 


1.001 


1.001 


1.007 


44 GHz 


















LFI24 


1.028 


1 


015 


1.007 


1.002 


1.000 


1.003 


1.009 


LFI25 


1.039 


1 


024 


1.013 


1.005 


1.000 


0.999 


1.000 


LFI26 


1.050 


1 


032 


1.017 


1.007 


1.000 


0.997 


0.997 


(C*) 44 


1.039 


1 


024 


1.012 


1.004 


1.000 


0.999 


1.002 


30 GHz 


















LFI27 


1.078 


1 


049 


1.026 


1.010 


1.000 


0.996 


0.998 


LFI28 


1.079 


1 


049 


1.026 


1.009 


1.000 


0.997 


1.002 


(C>30 


1.079 


1 


049 


1.026 


1.010 


1.000 


0.997 


1.000 



— SEVEM: Template fitting in map or wavelet domain; 

— Wi-fit: Template fitting in wavelet domain. 



Details of these methods may be found in Leach et al. 
(2008). These six algorithms make different assumptions 
about the data, and may use different combinations of fre- 
quency channels used as input. Comparing results from 
these methods (see Fig. 14) demonstrated the consistency 
of the CMB template and provided an estimate of the un- 
certainties in the reconstruction. A detailed comparison of 
the output of these methods, largely based on the CMB 
angular power spectrum, was used to select the CMB tem- 
plate that was removed from the frequency channel maps. 
The comparison was quantified using a jackknife procedure: 
each algorithm was applied to two additional sets of fre- 
quency maps made from the first half and second half of 
each pointing period. A residual map consisting of half the 
difference between the two reconstructed CMB maps was 
taken to be indicative of the noise level in the reconstruction 
from the full data set. The Needlet ILC (hereafter NILC) 
map was chosen as the CMB template because it had the 
lowest noise level at small scales. 

The CMB template was removed from the frequency 
channel maps after application of a filter in the spherical 
harmonic domain. The filter has a transfer function made of 
two factors. The first corresponds to the Gaussian beam of 
the channel to be cleaned; the second is a transfer function 
attenuating the multipoles of the CMB template that have 
low signal-to-noise ratio. It is designed in Wiener-like fash- 
ion, being close to unity up to multipoles around I = 1000, 
then dropping smoothly to zero with a cut-off frequency 
around i = 1700 (see Fig. 11). All angular frequencies above 
I = 3900 are completely suppressed. This procedure was 
adopted to avoid doing more harm than good to the small 
scales of the frequency channel maps where the signal-to- 
noisc ratio of the CMB is low. 



0.1 : 



0.01 r 




1000 2000 3000 

Fig. 11. Wiener- like filter function, plotted versus multi- 
pole, which was applied to produce the template for CMB 
removal. 

9.3. Description of Needlet ILC 

The NILC map was produced using the ILC method in the 
"needlet" domain. Needlets are spherical wavelets that al- 
low localisation both in multipole and sky direction. The 
input maps are decomposed into twelve overlapping mul- 
tipole domains (called "scales"), using the bandpass filters 
shown in Fig. 12 and further decomposed into regions of 
the sky. Independent ILCs are applied in each sky region 
at each needlet scale. Large regions are used at large scales, 
while smaller regions are used at fine scales. 

The NILC template was produced from all six HFI chan- 
nels, using the tight Galactic mask shown in figure 13, 
which covers 99.36% of the sky. Additional areas are ex- 
cluded on a per-channel basis to mask point sources. Future 
inclusion of the LFI channels will improve cleaning of low- 
frequency foregrounds such as synchrotron emission from 
the CMB template. Before applying NILC, pixels missing 
due to point source and Galactic masking are filled in by 
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3000 



Fig. 12. The bandpass filters, plotted versus multipolc, 
that define the spectral domains used in the NILC. 
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Fig. 14. Estimate of the rms error in the CMB subtraction. 
The map is histogram-equalised to bring out the details. 



Dispersion in the CMB rendition by different methods pro- 
vides an estimate of the uncertainties in the determination 
of the CMB, and thus in the subtraction process. The rms 
difference between the NILC map and the other CMB es- 
timates is shown in Fig. 14. As expected, the uncertainties 
are largest in the Galactic plane where the foregrounds to 
remove are strongest, and smallest around the Ecliptic poles 
where the noise levels are lowest. 



Fig. 13. Galactic mask used with NILC. 



a "diffusive inpainting" technique, which consists of replac- 
ing each missing pixel by the average of its neighbours and 
iterating to convergence. This is similar to solving the heat 
diffusion equation in the masked areas with boundary con- 
ditions given by the available pixel values at the borders of 
the mask. All maps are re-beamed to a common resolution 
of 5'. Re-beaming blows up the noise in the less resolved 
channels, but that effect is automatically taken into account 
by the ILC filter. 

The CMB template obtained after NILC processing is 
filtered to have the 'Wiener beam' shown in Fig. 11. The 
ILC coefficients are saved to be applied to the jackknifc 
maps for performance evaluation as described in Sect. 9.4.2 

9.4. Uncertainties in the CMB removal 

The uncertainties in the CMB removal have been gauged 
in two ways, firstly by comparing the CMB maps produced 
by the different algorithms and secondly by applying the 
NILC coefficients to jackknife maps. 

9.4.1. Dispersion of the CMB maps produced by the various 
algorithms. 

The methods that were used to produce the estimates 
of the CMB are diverse. They work by applying differ- 
ent algorithms (ILC, template fitting, or Bayesian pa- 
rameter estimation) in a variety of domains (pixel space, 
Needlet /wavelet space, or spherical harmonic coefficients). 
Each method carries out its optimisation in a different 
way and thus will respond to the foregrounds differently. 



9.4.2. CMB map uncertainties estimated by applying NILC 
filtering of jackknifes 

The cleanliness of the CMB template produced by the 
NILC filter can be estimated using jackknives. We apply 
the NILC filter to the maps built from the first and last 
halves of the ring set. The power distribution of the half- 
difference of the results provides us with a reliable estimate 
of the power of the noise in the NILC CMB template, (while 
previous results correspond to applying the NILC filter to 
the half-sum maps from which they can be derived). 

The jackknives allow estimates of the relative contri- 
butions of sky signal and noise to the total data power. 
Assume that the data are in the form X = S + N where 
S is the sky signal and N is the noise, independent of S. 
The total data power Var(X) decomposes as Var(AT) = 
Var(S') + Var(iV). One can obtain Var(iV) by applying the 
NILC filter to half difference maps, and Var(5) follows from 
Vai(X) — Var(iV). This procedure can be applied in pixel 
space, in harmonic space, or in pixel space after the maps 
have been bandpass-filtered, as described next. 

We first used pixel space jackknifing to estimate the 
spatial distribution of noise. Figure 15 shows a map of the 
local rms of the noise. We applied the NILC filter to a 
half-difference map and we display the square root of its 
smoothed squared values, effectively resulting in an esti- 
mate of the local noise rms. Using the same approach, we 
obtain an estimate of the angular spectrum of the noise 
in the NILC map, shown in Fig 16. That spectrum cor- 

1/2 

responds to an rms [(l/47r) J2 e (2£ + l)Cy of 11 [iK per 
pixel. The "features" in the shape of the noise angular spec- 
trum at large scale are a consequence of the needlet-based 
filtering (such features would not appear in a pixel-based 
ILC map). Recall that the coefficients of an ILC map are 
adjusted to minimize the total contamination by both fore- 
grounds and noise. The strength of foregrounds relative to 
noise being larger at coarse scales, the needlet-based ILC 
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the NILC CMB map. The colour scale is from to 30 fiK 
per pixel at resolution iV s id = 2048. 




Fig. 17. Local power of the NILC CMB template in the 
range £ = 500 ± 200. 

tends to let more noise in, with the benefit of better fore- 
ground rejection. 

The half-difference maps offer simple access to the 
power distribution of the residual noise in the estimated 
CMB template. However, it is more difficult to evaluate 
other residual contamination, since all fixed sky emissions 
cancel in half difference maps. Any such large-scale con- 
tamination is barely visible in the CMB template, since 
it is dominated by the CMB itself. However, contamina- 
tion is more conspicuous if one looks at intermediate scales. 
Figure 17 shows the local power of the CMB template af- 



Fig. 18. CMB-removed channel maps. From top to bottom, 
30, 44, and 70 GHz. The main galactic structures are clearly 
visible, as well as scanning strategy signatures at 44 and 
70 GHz. 

ter it is bandpassed to retain only multipolcs in the range 
I = 500 ± 200. This smooth version of the square of a band- 
passed map clearly shows where the errors in the compo- 
nent separation become large and so complicate some spe- 
cific science analyses. 

10. Infrastructure overview 

To organize the large number of data processing codes and 
data products, the DPC employs the Planck Integrated 
Data and Information System (IDIS). This allows flexi- 
ble development of the processing pipeline, while ensur- 
ing complete traceability and reproducibility of data prod- 
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ucts. For this, the most relevant components of IDIS 
are the Data Management Component (DMC) and the 
Process Coordinator (ProC), developed at the MPA Planck 
Analysis Centre (MPAC) at the Max-Planck-Institute for 
Astrophysics in Garching. Access to both components as 
well as to the Planck document and software manage- 
ment system is controlled via another IDIS component: 
the Federation Layer developed and maintained at the 
ESA ESTEC Research and Scientific Support Department 
(RSSD). 

Here we describe the essential features of the IDIS data 
processing components and their use at the DPC. A more 
detailed description of these components and their capabil- 
ities will be given in a future paper. 

10.1. Data Management Component 

The DMC organizes the storage and access to DPC 
data products. To combine optimal performance in data 
I/O with the data management capabilities of modern 
databases, scientific data are stored in files, while metadata 
identifying them are stored in a database. The data files 
can only be modified in synchronization with the database, 
preventing concurrent access to data objects via locking 
mechanisms. The DMC software supports several database 
management systems of various complexity; the LFI DPC 
operates an Oracle lOg database, which ensures good per- 
formance and stability. 

The DMC provides a uniform Application Programming 
Interface (API) for Fortran, C, C++, and Java, hiding all 
specific database operations from the user, who is therefore 
not required to have database experience. DMC data types 
are defined in the Data Definition Layer (DDL), which de- 
scribes data and metadata structures. The DDL supports 
inheritance of data types (e.g., a data type polar izedjnap 
can be inherited from a data type map) as well as associa- 
tion of data types (i.e., one data type containing a reference 
to another). 

In addition to the API, the DMC provides a Graphical 
User Interface (GUI), which supports user queries of the 
database and retrieval of information on the data. The 
GUI offers the user the ability to list the (meta-)data of 
specific objects and also to visualize the data in a simple 
way (although data can also be exported to other power- 
ful visualizing tools). The GUI also allows the display of 
history information on data objects, permitting the user to 
browse intermediate data products used in generating those 
objects, and the controlled deletion of data, observing de- 
pendencies of data types and maintaining the history infor- 
mation for all remaining data. For this, the DMC relies on 
additional metadata on the processing history of the data 
objects, which are generated by the Process Coordinator 
workflow engine (ProC). 

10.2. Pipeline management — the ProC workflow engine 

The ProC is a generic engine to construct, verify, and 
execute computational workflows. It comprises computing 
modules and data flows between them. The modules can be 
written in any programming language, provided they con- 
form to simple I/O format requirements described in an 
XML module description file. These interface files specify 
the input and output objects, as well as the parameters of 



the individual programs, in terms of DMC data types as 
described above. 

The ProC provides a Pipeline Editor to support graph- 
ical construction of data processing workflows. It allows 
users to arrange and connect computing modules of a work- 
flow in a clearly structured manner, and at the same time 
to configure the parameters of the algorithms used. It pro- 
vides control structures for data flow, for data object I/O 
and consistent parameter definition. 

The execution of workflows is controlled by a forward 
chaining algorithm, which ensures that modules are exe- 
cuted as soon as all necessary data products and parame- 
ters are known. If the same version of a module has been 
executed with identical inputs and parameters, the ProC 
will skip the execution and use the data product from the 
earlier execution for further processing. The ProC main- 
tains control of pipeline execution also on massive parallel 
computing environments. In the LFI implementation, the 
ProC communicates with the PBS (Portable Batch System) 
scheduling system to send jobs to the DPC cluster and to 
log their execution status. 

The ProC logs workflow executions on log files, which 
can also contain logging messages of the executed modules. 
Additionally, it creates so-called Pipeline-Run and Module- 
Run objects in the DMC, which are used to recover the 
generation history of data products (including versions of 
processing modules via MD5-sums). Besides the GUI, the 
ProC can also be executed from the command line. 

At the LFI DPC, the ProC is used to execute the official 
pipeline producing Planck data products. 

11. Discussion and conclusions 

We have described the status of the pipeline as it stands 
at the time of the ERCSC release and submission of the 
Planck early papers. All the algorithms run during this pro- 
cess have been verified, validated, and tested before launch 
and the start of operations using realistic simulations. This 
allowed us to begin analyzing the data as soon as they were 
acquired from the first day of operations. The entire Level 1 
pipeline suffered no significant problems, and all of the data 
were transformed efficiently from telemetry packets to time- 
lines. At present, the Level 2 pipeline is capable of provid- 
ing relative calibration to an overall statistical accuracy in 
the range 0.05-0.1% and absolute calibration at around the 
1% level. The beams are accurately characterised down to 
— lOdB. We expect to improve many aspects in the near fu- 
ture. Concerning the calibration, our intention is to reach 
the levels determined by the stability of the instrument. For 
the beam reconstruction, our aim is to improve the char- 
acterisation of the far side lobes and to refine the entire 
beam reconstruction pipeline, with particular attention to 
polarization measurements. 
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